1 Méthodologie

Les traitements statistiques ont été réalisés avec le logiciel R (version 4.1.1). Par ailleurs, les opérations géostatistiques ont été réalisées à l’aide du package ‘gstat’ disponible sur R. Les transformations Box-Cox sont réalisée à l’aide de la librairie ‘AID’. Les graphes ont quant à été générés avec l’extension ‘ggplot’.

2 Traitement des données

Les données utilisées ont été récoltées par plusieurs laboratoires en Wallonie. Le set que nous utilisons est un extrait de ce large jeu de données comprenant les concentrations trois éléments traces métalliques (ETM) pour une province Wallonne.
Plus précisément, il s’agit de concentrations en Plomb, Zinc et Cuivre récoltés à de multiples sites d’échantillonnage contenue dans la province de Liège.

mydata <- fread('data.csv')
summary(mydata)
##      FID_1             X                Y               ID           
##  Min.   :  290   Min.   :195331   Min.   : 93704   Length:2541       
##  1st Qu.: 8164   1st Qu.:217492   1st Qu.:127264   Class :character  
##  Median :13670   Median :224149   Median :138824   Mode  :character  
##  Mean   :12669   Mean   :229780   Mean   :137604                     
##  3rd Qu.:16992   3rd Qu.:240142   3rd Qu.:147500                     
##  Max.   :22990   Max.   :291500   Max.   :163250                     
##                                                                      
##        Pb                 Zn                 Cu            SIGLE_PEDO       
##  Min.   :   8.127   Min.   :    5.80   Min.   :   2.016   Length:2541       
##  1st Qu.:  26.135   1st Qu.:   85.12   1st Qu.:  12.905   Class :character  
##  Median :  32.000   Median :  112.53   Median :  15.485   Mode  :character  
##  Mean   :  52.612   Mean   :  270.92   Mean   :  22.576                     
##  3rd Qu.:  48.916   3rd Qu.:  182.60   3rd Qu.:  20.000                     
##  Max.   :1149.970   Max.   :25364.31   Max.   :2234.267                     
##  NA's   :1038       NA's   :157        NA's   :273                          
##    region_etm           GS       
##  Min.   : 2.000   Min.   : 3.00  
##  1st Qu.: 6.000   1st Qu.: 9.00  
##  Median : 7.000   Median : 9.00  
##  Mean   : 6.749   Mean   :12.57  
##  3rd Qu.: 8.000   3rd Qu.:16.00  
##  Max.   :17.000   Max.   :43.00  
## 
head(mydata)
##    FID_1      X      Y         ID       Pb        Zn       Cu SIGLE_PEDO
## 1:   409 238401 153664 02/LA00951 32.14414  10.97368 13.85714       Aba1
## 2:   654 237838 152974 03/LA01232 18.59821        NA 12.00000       Aba1
## 3:   511 224562 137188 02/LA04725 32.14414 135.18421 14.80952       Aba1
## 4:   473 214050 132070 02/LA04323 40.25225  57.55263 12.90476       Aca0
## 5:   544 214043 133536 02/LA05640       NA 175.97368 10.04762       Ada0
## 6:   512 221460 135511 02/LA04727 45.65766 252.28947 15.76190       Ada1
##    region_etm GS
## 1:          6  9
## 2:          6  9
## 3:          6  9
## 4:          6  9
## 5:          6  9
## 6:          6  9

Le jeu de données contient 10 variables à savoir : un identifiant, des coordonnées (latitude et longitude), 3 éléments traces métalliques et 3 autres variables dont SIGLE-PEDO qui informe sur le type de sol.

Description du jeu de données
Description du jeu de données

2.1 Carte de la position des points de mesures des ETM et limite de la province étudiée.

# Create grid (using a margin to make sure you englobe the whole province)
gridsize = 1000
margin = 5000
x <- seq(floor(min(mydata$X)-margin), # from minimum longitude
         ceiling(max(mydata$X+margin)), # to maximum longitude
         by=gridsize)
y <- seq(floor(min(mydata$Y)-margin), # from minimum latitude
         ceiling(max(mydata$Y)+margin), # to maximum latitude
         by=gridsize)
grid <- as.data.table(expand.grid(X=x, Y=y)) 

# Create a spatial version of your grid
mydataSpatial <- copy(grid)
coordinates(mydataSpatial) <- ~X+Y
proj4string(mydataSpatial) <- CRS("+init=epsg:31370") # Specify coordinate system (Lambert belge 1972)
# Load provinces shapefile and specify it's coordinate system
prov <- readOGR('BEL_ADM2.shp', use_iconv = TRUE, encoding = "UTF-8")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\nicd\Desktop\Portfolio\templates\project\mod_spatiale\BEL_ADM2.shp", layer: "BELGIUM_-_Provinces"
## with 11 features
## It has 8 fields
prov <- spTransform(prov, CRS("+init=epsg:31370"))
plot(prov) # Display provinces

prov$NAME_2
##  [1] "Antwerpen"       "Bruxelles"       "Vlaams Brabant"  "Brabant Wallon" 
##  [5] "West-Vlaanderen" "Oost-Vlaanderen" "Hainaut"         "Liège"          
##  [9] "Limburg"         "Luxembourg"      "Namur"
# Transform grid coordinates system to match provinces
mydataSpatial <- spTransform(mydataSpatial, CRS("+init=epsg:31370"))
## Warning: PROJ support is provided by the sf and terra packages among others
# Get index of points in the province of Liège
prov.id <- over(prov[prov$NAME_2 == "Liège", ],
               mydataSpatial, 
               returnList = TRUE)

# Select rows in your original grid that are located in Liège
prov.grid = grid[prov.id[[1]],]
# Display selected grid points in red and data location
ggplot()+
  geom_point(data=prov.grid, aes(x=X, y=Y), color='black', shape = 3)+
  geom_point(data=grid, aes(x=X, y=Y), color='black', shape = 3) +
 geom_point(data=prov.grid, aes(x=X, y=Y), color='red', shape = 3)+
 geom_point(data=mydata,aes(x=X,y=Y,shape="Point de mesure"),color="black")+
 ggtitle("Limite de la province étudiée (Liège) et positions des points de mesures")+
 xlab("Longitude (Lambert belge 1972)")+
 ylab("Latitude (Lambert belge 1972)")+
 scale_shape_manual("",values=15)

2.2 Nettoyage et visualisation des échantillons par ETM

2.2.1 Pb sample

Pb <- na.omit(subset(mydata,select=c("X","Y","Pb")))

ggplot() +
  geom_point(data=Pb, aes(x=X, y=Y, color=Pb), size=2) +
  scale_color_gradient(low='blue',high ='yellow')+
  xlab("Longitude") +
  ylab("Latitude") +
  ggtitle("Sample Pb")

2.2.2 Zn sample

Zn <- na.omit(subset(mydata,select=c("X","Y","Zn")))

#Plot

ggplot() +
  geom_point(data=Zn, aes(x=X, y=Y, color=Zn), size=2) +
  scale_color_gradient(low='blue',high ='yellow')+
  xlab("Longitude") +
  ylab("Latitude") +
  ggtitle("Sample Zn")

2.2.3 Cu sample

Cu <- na.omit(subset(mydata,select=c("X","Y","Cu")))

#Plot

ggplot() +
  geom_point(data=Cu, aes(x=X, y=Y, color=Cu), size=2) +
  scale_color_gradient(low='blue',high ='yellow')+
  xlab("Longitude") +
  ylab("Latitude") +
  ggtitle("Sample Cu")

2.3 Analyse exploratoire des données

Le jeu de données des ETM ayant été l’objet d’un premier traitement, les statistiques de bases ont été calculées.

mydata %>%
  select(Pb, Zn, Cu) %>%
  summary %>%
  pander()
Pb Zn Cu
Min. : 8.127 Min. : 5.80 Min. : 2.016
1st Qu.: 26.135 1st Qu.: 85.12 1st Qu.: 12.905
Median : 32.000 Median : 112.53 Median : 15.485
Mean : 52.612 Mean : 270.92 Mean : 22.576
3rd Qu.: 48.916 3rd Qu.: 182.60 3rd Qu.: 20.000
Max. :1149.970 Max. :25364.31 Max. :2234.267
NA’s :1038 NA’s :157 NA’s :273 
chart.Correlation(mydata[,.(Zn,Pb,Cu)]) 
## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

Etant donné que les distributions observées des ETM ne s’apparente pas à des lois normales, il est nécessaire d’appliquer une transformation à chacune des trois distributions pour pouvoir utilisé les modèles.

Une forte corrélation, soit la meilleure, entre le plomb et le zinc est observée, avec un coefficient de l’ordre de 0.73. Le Zinc est presqu’aussi bien corrélé au cuivre avec un coefficient de 0.73. Enfin, la corrélation entre le plomb et le cuivre est la plus faible, soit un coefficient de 0.59

2.4 Correction de la distribution

Il est nécessaire, pour satisfaire aux hypothèses des modèles uitilisé pour les analyses et prédictions d’avoir des données suivant une distribution normale. La transformation Box-Cox a été utilisée afin de corriger la distribution des ETM lorsque celle-ci ne fut pas similaire à une distribution normale.

Pour une variable Z à valeurs non négatives, la transformée de Box-Cox est définit de la manière suivante :

\[Y_i = \begin{cases} \frac{Z^{\theta}_{i}-1}{\theta}, & \text{si } \theta \neq 0 \\ ln(Z_{i}), & \text{si } \theta = 0 \end{cases}\]

Le paramètre θ est alors choisi pour que la variable transformée Yi s’approche au mieux d’une distribution normale.

2.4.1 Box-cox sur la distribution des échantillons de Plomb

# Plomb
BPb <- boxcoxnc(Pb$Pb, verbose = FALSE)

Pb$Pbbx <- BPb$tf.data

2.4.2 Box-cox sur la distribution des échantillons de Zinc

BZn <- boxcoxnc(Zn$Zn, verbose = FALSE)

Zn$Znbx <- BZn$tf.data

2.4.3 Box-cox sur la distribution des échantillons de Cuivre

BCu <- boxcoxnc(Cu$Cu, verbose = FALSE)

Cu$Cubx <- BCu$tf.data

3 Analyse et modélisation de la dépendence spatiale

Afin de mesurer et modéliser la dépendance spatiale qu’il pourrait y avoir entre les valeurs d’un ETM, la fonction espérance et le variogramme sont utilisés. La fonction espérance permet dans un premier temps de rendre compte de la présence ou non de stationnarité d’ordre 1. Dans la négative, l’espérance est donc fonction des coordonnées et le paramètre β alors être estimé à l’aide d’un modèle de régression linéaire (OLS, WLS et GLS). Le variogramme estime la dépendance dans l’espace en mesurant la dissimilitude des données en fonction de la distance h, avec :

\[ \gamma(h)= \frac{1}{2} Var[Z(x+h) - Z(x)] \quad \forall x \in D \] Où - Z est la variable spatiale étudiée - x est un vecteur de coordonnées - h est le vecteur de distance - D est le domaine géographique considéré

On part de l’hypothèse (courramment faite) : hmax <

La figure ci-dessous illustre le variogramme et les paramètres qu’il comporte.

variogramme
variogramme

3.1 Espérances et variogramme Plomb (Pb)

3.1.1 Espérence du Plomb B-C

Pbbx.lm <- lm(Pbbx~X+Y,data=Pb)
scatter3d(x=Pb$X,z=Pb$Y,y=Pb$Pbbx,
          xlab = 'Longitude', zlab = 'Lattitude', ylab = 'Pbbx')
## Le chargement a nécessité le package : mgcv
rglwidget(width = 300, height = 300)

Non stationnaire d’ordre 1 + on voit des outliers

3.1.2 Suppression de la tendance et outlier

# Suppression tendance
Pb$Pbbx_res <- Pb$Pbbx-(predict(Pbbx.lm,Pb))
scatter3d(x=Pb$X,z=Pb$Y,y=Pb$Pbbx_res,
          xlab='Longitude',zlab='Latitude',ylab='Pbbx_res')
rglwidget(width = 300,height = 300)
# Détection potentiel Outlier
res.Sd <- sd(Pb$Pbbx_res)
res.Mean <- mean(Pb$Pbbx_res)
out <- nrow(Pb) - nrow(Pb[Pbbx_res %between% c(res.Mean-3*res.Sd,res.Mean+3*res.Sd),])
cat("nombre d'outlier potentiel =", out)
## nombre d'outlier potentiel = 15
# Suppression outlier
Pb <- Pb[Pbbx_res %between% c(res.Mean-3*res.Sd,res.Mean+3*res.Sd),]

3.1.3 Variogramme résidu du Plomb Bx

Pb.gstat <- gstat(formula=Pbbx_res~1, data=Pb,locations = ~X+Y)
Pb.vario <- variogram(Pb.gstat,cutoff=40000,width=3000)
vg.exp <- vgm(model = 'Sph',nugget = 2e-4)
fit.Pb <- fit.variogram(Pb.vario, model = vg.exp,fit.method = 6)
plot(Pb.vario,model=fit.Pb, main="Variogram résidus du Plomb Bx",pch=16,col='black')
trellis.focus('panel',1,1)
llines(x=c(0,max(Pb.vario$dist)),y=c(var(Pb$Pbbx_res)),col = 'red',lty=2)
## NULL
trellis.unfocus()

pander(fit.Pb)
model psill range kappa ang1 ang2 ang3 anis1 anis2
Nug 0.0001076 0 0 0 0 0 1 1
Sph 0.0004151 20021 0.5 0 0 0 1 1

3.2 Espérances et variogramme Zinc (Zn)

3.2.1 Espérence du Zinc

Znbx.lm <- lm(Znbx~X+Y,data=Zn)
scatter3d(x=Zn$X,z=Zn$Y,y=Zn$Znbx,
          xlab = 'Longitude', zlab = 'Lattitude', ylab = 'Znbx')
rglwidget(width = 300, height = 300)

Non stationnaire d’ordre 1 + on voit des outliers

3.2.2 Suppression de la tendance et outlier

# Suppressions tendance
Zn$Znbx_res <- Zn$Znbx-(predict(Znbx.lm,Zn))
scatter3d(x=Zn$X,z=Zn$Y,y=Zn$Znbx_res,
          xlab='Longitude',zlab='Latitude',ylab='Znbx_res')
rglwidget(width = 300,height = 300)
# Détection potentiel Outlier
res.Sd <- sd(Zn$Znbx_res)
res.Mean <- mean(Zn$Znbx_res)
out <- nrow(Zn) - nrow(Zn[Znbx_res %between% c(res.Mean-3*res.Sd,res.Mean+3*res.Sd),])
cat("nombre d'outlier potentiel =", out)
## nombre d'outlier potentiel = 47
# Suppression outlier
Zn <- Zn[Znbx_res %between% c(res.Mean-3*res.Sd,res.Mean+3*res.Sd),]

3.2.3 Variogramme résidu du Zinc Bx

Zn.gstat <- gstat(formula=Znbx_res~1, data=Zn,locations = ~X+Y)
Zn.vario <- variogram(Zn.gstat,cutoff=40000,width=3000)

vg.exp <- vgm(model = 'Sph',nugget = 0.002)
fit.Zn <- fit.variogram(Zn.vario, model = vg.exp,fit.method = 6)
plot(Zn.vario,model=fit.Zn ,main="Variogram résidus du Zinc",pch=16,col='black')

trellis.focus('panel',1,1)
llines(x=c(0,max(Zn.vario$dist)),y=c(var(Zn$Znbx_res)),col = 'red',lty = 2)
## NULL
trellis.unfocus()

pander(fit.Zn)
model psill range kappa ang1 ang2 ang3 anis1 anis2
Nug 0.001452 0 0 0 0 0 1 1
Sph 0.007236 17386 0.5 0 0 0 1 1

3.3 Espérances et variogramme Cuivre (Cu)

3.3.1 Espérence du Cuivre

Cubx.lm <- lm(Cubx~X+Y,data=Cu)
scatter3d(x=Cu$X,z=Cu$Y,y=Cu$Cubx,
          xlab = 'Longitude', zlab = 'Lattitude', ylab = 'Cubx')
rglwidget(width = 300, height = 300)

Non stationnaire d’ordre 1 + on voit des outliers

3.3.2 Suppression de la tendance et outlier

# Suppressions tendance
Cu$Cubx_res <- Cu$Cubx-(predict(Cubx.lm,Cu))
scatter3d(x=Cu$X,z=Cu$Y,y=Cu$Cubx_res,
          xlab='Longitude',zlab='Latitude',ylab='Cubx')
rglwidget(width = 300,height = 300)
# Détection potentiel Outlier
res.Sd <- sd(Cu$Cubx_res)
res.Mean <- mean(Cu$Cubx_res)
out <- nrow(Cu) - nrow(Cu[Cubx_res %between% c(res.Mean-3*res.Sd,res.Mean+3*res.Sd),])
cat("nombre d'outlier potentiel =", out)
## nombre d'outlier potentiel = 34
# Suppression outlier
Cu <- Cu[Cubx_res %between% c(res.Mean-3*res.Sd,res.Mean+3*res.Sd),]

3.3.3 Variogramme résidu du Cuivre Bx

Cu.gstat <- gstat(formula=Cubx_res~1, data=Cu,locations = ~X+Y)
Cu.vario <- variogram(Cu.gstat,cutoff=40000,width=3000)

vg.exp <- vgm(model = 'Sph',nugget = 0.01)
fit.Cu <- fit.variogram(Cu.vario, model = vg.exp,fit.method = 6)
plot(Cu.vario,model=fit.Cu, main="Variogram résidus du Cuivre",pch=16,col='black')
trellis.focus('panel',1,1)
llines(x=c(0,max(Cu.vario$dist)),y=c(var(Cu$Cubx_res)),col = 'red',lty=2)
## NULL
trellis.unfocus()

pander(fit.Cu)
model psill range kappa ang1 ang2 ang3 anis1 anis2
Nug 0.01108 0 0 0 0 0 1 1
Sph 0.01495 15179 0.5 0 0 0 1 1

4 Prédiction des ETM

Demande normalité, stationnarité d’ordre 1 et 2

plot.predictions <- function(pred,data, plot.title,legend){
  ggplot() + 
  geom_tile(data=pred, 
            aes(x = X, y = Y, fill = predi)) +
  geom_point(data=data, 
             aes(x=X, y=Y, color="Measurement points"),
             shape=18,
             size=2) +
  scale_color_manual("", values="black") +
  scale_fill_gradientn(name=legend, colors=c('royalblue','green3', 'yellow', 'red')) +
  theme(legend.key = element_rect(fill = "green3", 
                                  color = NA)) +
  xlab("Longitude") +
  ylab("Latitude") +
  ggtitle(plot.title)
}

4.1 Prédiction par la méthode de distance inverse

La prédiction par distance inverse (Inverse distance weighting [IDW]) est un prédicteur exact qui possède l’avantage de prendre en compte la distance entre le point prédit et les échantillons. On attribue ainsi à chaque points xi poids qui est proportionnel à l’inverse de la distance entre ce point et x0 ou l’on réalise la prédiction. Formellement, ce poids est donné par :

\[\lambda_{i} = \frac{\frac{1}{||h_{0i}||^{\theta}}}{\sum^{n}_{j=1}\frac{1}{||h_{0j}||^{\theta}}} \quad avec : \theta > 0\] Le paramètre θ détermine donc la vitesse de diminution du poids d’un point en fonction de la distance ||h||. Ce paramètre est optimisé en minimisant la somme du carré des écarts lors de la leave-one-out cross validation (LOOCV).

4.1.1 Plomb

4.1.1.1 Optimisation de θ par LOOCV

powers <- seq(0.5,5,0.5)
pmse <- data.table(power = powers, mse = rep(0,length(powers)) )# Power mean squared error
for (p in powers){
  pse <- rep(0,nrow(Pb)) # Power squared errors
  for (i in 1:nrow(Pb)){
     point.idw <- idw(formula = Pbbx~1,
              data = Pb[-i,],
              locations = ~X+Y,
              newdata = Pb[i,],
              idp = p,
              nmax=20,
              maxdist=30000,
              debug.level = 0)
     pse[i] <- (point.idw$var1.pred - Pb$Pbbx[i])^2
  }
  pmse[power==p,"mse"] = mean(pse)
}

plot(pmse)

4.1.1.2 Plot résultat IDW Plomb

#Plot résultat
Pb.idw <- data.table(idw(formula = Pbbx~1,data=Pb,location=~X+Y,newdata=prov.grid,idp=1.5,nmax=20))
## [inverse distance weighted interpolation]
setnames(Pb.idw,"var1.pred","predBx")

#Transfo inverse
alpha = BPb$lambda.hat
Pb.idw[,predi := exp(log(alpha * predBx + 1) / alpha)]

plot.predictions(Pb.idw,Pb,"Concentration en plomb - Prediction Distance inverse","Pb predictions")

4.1.2 Zinc

4.1.2.1 Optimisation de θ par LOOCV

powers <- seq(0.5,5,0.5)
pmse <- data.table(power = powers, mse = rep(0,length(powers)) )# Power mean squared error
for (p in powers){
  pse <- rep(0,nrow(Zn)) # Power squared errors
  for (i in 1:nrow(Zn)){
     point.idw <- idw(formula = Znbx~1,
              data = Zn[-i,],
              locations = ~X+Y,
              newdata = Zn[i,],
              idp = p,
              nmax=20,
              maxdist=30000,
              debug.level = 0)
     pse[i] <- (point.idw$var1.pred - Zn$Znbx[i])^2
  }
  pmse[power==p,"mse"] = mean(pse)
}

plot(pmse)

4.1.2.2 Plot résultat IDW Zinc

Zn.idw <- data.table(idw(formula = Znbx~1,data=Zn,location=~X+Y,newdata=prov.grid,idp=1,nmax=20,maxdist=30000))
## [inverse distance weighted interpolation]
setnames(Zn.idw,"var1.pred","predBx")

# Transformation inverse
alpha = BZn$lambda.hat
Zn.idw[,predi := exp(log(alpha * predBx + 1) / alpha)]

plot.predictions(Zn.idw,Zn,"Concentration en Zinc - Prediction Distance inverse","Zn predictions")

4.1.3 Cuivre

4.1.3.1 Optimisation de θ par LOOCV

powers <- seq(0.4,2,0.2)
pmse <- data.table(power = powers, mse = rep(0,length(powers)) )# Power mean squared error
for (p in powers){
  pse <- rep(0,nrow(Cu)) # Power squared errors
  for (i in 1:nrow(Cu)){
     point.idw <- idw(formula = Cubx~1,
              data = Cu[-i,],
              locations = ~X+Y,
              newdata = Cu[i,],
              idp = p,
              nmax=20,
              maxdist=30000,
              debug.level = 0)
     pse[i] <- (point.idw$var1.pred - Cu$Cubx[i])^2
  }
  pmse[power==p,"mse"] = mean(pse)
}

plot(pmse)

4.1.3.2 Plot résultat IDW Cuivre

Cu.idw <- data.table(idw(formula = Cubx~1,data=Cu,location=~X+Y,newdata=prov.grid,idp=1.25,nmax=20,maxdist=30000))
## [inverse distance weighted interpolation]
setnames(Cu.idw,"var1.pred","predBx")

# Transformation inverse
alpha = BCu$lambda.hat
Cu.idw[,predi := exp(log(alpha * predBx + 1) / alpha)]

plot.predictions(Cu.idw,Cu,"Concentration en Cuivre - Prediction Distance inverse","Cu predictions")

4.2 Via le krigeage

Le Krigeage ou BLUP (Best Linear Unbiased Predictor) est une technique géostatistique de modélisation spatiale qui présente l’avantage de prendre en compte les distances entre les données et le point estimé et la structure spatiale (grâce au variogramme). Ainsi, le krigeage se base sur une hypothèse d’autocorrélation spatiale des données en posant que deux données proches dans l’espace possèdent des caractéristiques similaires. Le krigeage a basiquement un modèle similaire à la régression linéaire mais les erreurs y sont supposées dépendantes spatialement :

\[Z(x) = \mu(x) + \delta(x) \quad \forall x \in D\] Où : - est l’espérance de Z - est une fonction aléatoire stationnaire

4.2.1 Plomb

doublons <- which(duplicated(Pb$X,Pb$Y))
Pb2<-Pb[-doublons,]
Pb.krig <- data.table(krige(formula = Pbbx_res~1, data = Pb2,locations = ~X+Y,
                            newdata = prov.grid, model = fit.Pb))
## [using ordinary kriging]
setnames(Pb.krig, c("var1.pred", "var1.var"), c("pred", "Pb.varkrig"))
Pb.krig$pred <- Pb.krig$pred+predict(Pbbx.lm,Pb.krig)
alpha <- BPb$lambda.hat
Pb.krig[,predi := exp(log(alpha * pred + 1) / alpha)]

plot.predictions(Pb.krig,Pb,"Concentration en Plomb - Krigeage","Pb prediction")

4.2.2 Zinc

doublons <- which(duplicated(Zn$X,Zn$Y))
Zn2<-Zn[-doublons,]
Zn.krig <- data.table(krige(formula = Znbx_res~1, data = Zn2,locations = ~X+Y,
                            newdata = prov.grid, model = fit.Zn))
## [using ordinary kriging]
setnames(Zn.krig, c("var1.pred", "var1.var"), c("pred", "Zn.varkrig"))
Zn.krig$pred <- Zn.krig$pred+predict(Znbx.lm,Zn.krig)
alpha <- BZn$lambda.hat
Zn.krig[,predi := exp(log(alpha * pred + 1) / alpha)]
plot.predictions(Zn.krig,Zn,"Concentration en Zinc - Krigeage","Zn predictions")

4.2.3 Cuivre

doublons <- which(duplicated(Cu$X,Cu$Y))
Cu2<-Cu[-doublons,]

Cu.krig <- data.table(krige(formula = Cubx_res~1, data = Cu2,locations = ~X+Y,
                            newdata = prov.grid, model = fit.Cu,nmax=1))
## [using ordinary kriging]
setnames(Cu.krig, c("var1.pred", "var1.var"), c("pred", "Cu.varkrig"))
Cu.krig$pred <- Cu.krig$pred+predict(Cubx.lm,Cu.krig)
alpha <- BCu$lambda.hat
Cu.krig[,predi := exp(log(alpha * pred + 1) / alpha)]
plot.predictions(Cu.krig,Cu,"Concentration en Cuivre - Krigeage","Cu predictions")

4.3 Via cokrigeage de Pb avec Zn

Dans le cas ou deux champs aléatoires sont corrélés, le krigeage peut être précisé en combinant ces deux échantillons. La démarche est alors similaire à celle du Krigeage défini ci-dessus mais prend alors en compte une seconde variable pour la prédiction désirée. L’amélioration de la prédiction est d’autant plus importante lorsque : - Les variables sont fortement corrélées ; - La taille de la seconde variable est grande (par rapport à celle qu’on veut prédire) ; - Les positions échantillonnées de la seconde variable varient de ceux de la première.

g <- gstat(id="Pb", formula=Pbbx~1, data=Pb2, locations= ~X+Y)
g <- gstat(g,id='Zn',formula = Znbx~1, data = Zn2, locations = ~X+Y)
v.cross <- variogram(g, cutoff=40000, width=3000)
plot(v.cross)

LMC <- fit.lmc(v.cross,g,model = fit.Pb, correct.diagonal = 1.1,fit.method=6)
LMC
## data:
## Pb : formula = Pbbx`~`1 ; data dim = 1077 x 3
## Zn : formula = Znbx`~`1 ; data dim = 1530 x 3
## variograms:
##          model       psill    range
## Pb[1]      Nug 0.001364857     0.00
## Pb[2]      Sph 0.001252742 20020.83
## Zn[1]      Nug 0.002195824     0.00
## Zn[2]      Sph 0.011670383 20020.83
## Pb.Zn[1]   Nug 0.001573801     0.00
## Pb.Zn[2]   Sph 0.003476007 20020.83
## ~X + Y
g <- copy(LMC)
g <- gstat(g, id='Zn', formula = Znbx~1,data = Zn2,locations = ~X+Y,model = LMC$model$Zn)
plot(v.cross,model=g$model,col='black',pch=16)

Pb.cok <- data.table(predict(g,prov.grid))
## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]
setnames(Pb.cok, c('Pb.pred','Pb.var'),c('pred','Pb.varcok'))
#Pb.cok$pred <- Pb.cok$pred+predict(Pbbx.lm,Pb.cok)
alpha <- BPb$lambda.hat
Pb.cok[,predi := exp(log(alpha * pred + 1) / alpha)]
plot.predictions(Pb.cok,Pb,'Concentration en Plomb - Cokrigeage','Prediction')

4.4 Simulation conditionnelles du plomb

Dans le cas d’une simulation conditionnelle, on dispose et prend en compte les valeurs connues des échantillons. Les simulations sont alors basées sur le variogramme γ(||h||) et conditionnées par les concentrations connues par échantillonnages. En réalisant plusieurs de ces simulations pour le Plomb, il sera possible de définir les zones possiblement contaminée dont la concentration dépasse un seuil fixé. Les zones qui dépassent ce seuil dans plus de x pourcent des simulations peuvent alors être suspectées d’être contaminées.